
********************************************************************************
********************************************************************************
********************************************************************************
********* AVERAGE TREATMENT EFFECTS ON THE TREATED (ATT) ***********************

clear all
do "$maindir/Code/bootstrapres.do"
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"


********************************************************************************
************************ REGRESSIONS *******************************************
********************************************************************************

** regressions for average treatment effects on treated - standard econometrics
reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=., absorb(i.Household i.month_year) vce(cluster Household)
nlcom exp(_b[1.treated])-1, post
est sto att_simple
gen newsamp = e(sample) /* to keep stable sample across regressions */

reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=. & newsamp==1 , absorb(i.Household#i.end_month i.month_year) vce(cluster Household)
nlcom exp(_b[1.treated])-1, post
est sto att_hhmonth

reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=. & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
nlcom exp(_b[1.treated])-1, post
est sto att_hhmonthcounty


* regressions in levels
reghdfe total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household i.month_year) vce(cluster Household)
est sto att_simple_levels

reghdfe total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year) vce(cluster Household)
est sto att_hhmonth_levels

reghdfe total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
est sto att_hhmonthcounty_levels


* checking parallel trends 
reghdfe log_total_mmbtu ib12.months_since_treat c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)

matrix coefs = (.)
matrix vars = (.)
forval i = 1/24 {
	nlcom exp(_b[`i'.months_since_treat])-1
	matrix coef = r(b)
	matrix var = r(V)
	scalar nobs = e(N)
	matrix colnames coef = "`i'.months_since_treat"
	matrix colnames var = "`i'.months_since_treat"
	
	matrix coefs = (coefs, coef)
	matrix vars = (vars, var)
}

matrix coefs = coefs[1..., 2...]
matrix vars = vars[1..., 2...]
matrix vars = diag(vars)

bootstrapres coefs vars nobs
est sto att_eventstudy

coefplot (att_eventstudy), keep(*months_since_treat) baselevels  ///
			vertical label ytitle("Estimated Energy Savings") ///
			xtitle("Months From Treatment (Weatherization)") xline(12.5) ///
			xlabel(, labsize(small)) yline(0) ///
			graphregion(color(white)) bgcolor(white)
graph export "$maindir/Results/Figures/ATT_paralleltrends.png", replace


****** separating effects for gas versus electricity usage
reghdfe log_gas_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
nlcom exp(_b[1.treated])-1, post
est sto att_loggas

reghdfe log_electric_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
nlcom exp(_b[1.treated])-1, post
est sto att_logelectric

reghdfe gas_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
est sto att_gaslevels

reghdfe electric_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) vce(cluster Household)
est sto att_electriclevels
 
esttab att_hhmonthcounty att_loggas att_logelectric att_hhmonthcounty_levels att_gaslevels att_electriclevels ///
				using "$maindir/Results/Tables/ATT_gasVSelec.tex", ///
				replace b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)  ///
				rename(_nl_1 a 1.treated a) ///
				keep(a) coeflabels(a "IHWAP Treatment") nonotes ///
				stats(N, fmt(%9.0gc) labels("Observations"))



************************** ML average treatment effects
*** ml pct savings
egen totuse = total(pred_full) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & tau_ml!=. & newsamp==1 
egen totsave = total(tau_ml) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & tau_ml!=. & newsamp==1 
gen pct_save = totsave/totuse

sum pct_save
matrix ml_ate = (`r(mean)')
sca nobs = `r(N)'
gen mlsamp = newsamp if pct_save!=.
replace mlsamp = 0 if mlsamp==.

* bootstrapped SEs for ml pct savings
matrix input ml_ates = (.)
forval i = 1/200 {
    di "Starting iteration `i'"
	qui: gen weightuse`i' = pred_boots`i'*weight_boots`i'
	qui: egen totuse`i' = total(weightuse`i') if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & tau`i'!=. & newsamp==1 
	qui: gen weightsave`i' = tau`i'*weight_boots`i'
	qui: egen totsave`i' = total(weightsave`i') if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & tau`i'!=. & newsamp==1 
	qui: gen pct_save`i' = totsave`i'/totuse`i'
	qui : sum pct_save`i'
	qui: sca avgsaved`i' = `r(mean)'
	qui: matrix ml_ates = (ml_ates \ avgsaved`i')
	qui: drop weightuse`i' weightsave`i' totuse`i' totsave`i' pct_save`i'
}
mata : st_matrix("ml_ate_sd", sqrt(variance(st_matrix("ml_ates"))))
mata : st_matrix("ml_ate_var", variance(st_matrix("ml_ates")))

prog def bootstraplab, eclass
	args b V N samp /* input coef matrix, variance matrix, and number of obs */
	mat colnames `b'="treated"
	mat colnames `V'="treated"
	mat rownames `V'="treated"
	scalar obs = `N'
	ereturn post `b' `V', obs(`=scalar(obs)') esample(`samp') 
	ereturn local cmd "bootstraplab"
	ereturn local properties "b V"
	ereturn display
end
bootstraplab ml_ate ml_ate_var nobs
est store ml_simple


* ml savings in levels
sum tau_ml if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & newsamp==1 
matrix ml_ate = (`r(mean)')
sca nobs = `r(N)'

* bootstrapped SEs for ml savings in levels
matrix input ml_ates = (.)
forval i = 1/200 {
    di "Starting iteration `i'"
	qui : sum tau`i' [fweight = weight_boots`i'] if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & newsamp==1 
	sca avgsaved`i' = `r(mean)'	
	matrix ml_ates = (ml_ates \ avgsaved`i')
}
mata : st_matrix("ml_ate_sd", sqrt(variance(st_matrix("ml_ates"))))
mata : st_matrix("ml_ate_var", variance(st_matrix("ml_ates")))
bootstraplab ml_ate ml_ate_var nobs
est store ml_simple_levels


**** "Event study" graphs for ML method
**** comparing true usage versus predicted usage
reg total_mmbtu ibn.months_since_treat if twoyears_prepost==1 & meterdays_monthly!=. & tau_ml!=.  & newsamp==1  , vce(cluster Household) nocons
est sto totaluse

*** will use cross-validated usage for pre-treatment
gen pred_cv = pred_full
replace pred_cv = cvpreds if treated==0 & pred_full!=.

reg pred_cv ibn.months_since_treat if twoyears_prepost==1 & meterdays_monthly!=. & tau_ml!=.  & newsamp==1 , nocons

* Bootsatrpping standard errors for event study graph
matrix ml_ate_mst = e(b)
scalar nobs = e(N)
matrix input ml_mst = (., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., ., .,)
forval i = 1/200 {
	qui: reg pred_boots`i' ibn.months_since_treat [fweight = weight_boots`i'] if twoyears_prepost==1 & meterdays_monthly!=. & tau`i'!=., nocons
	matrix tempmat = e(b)
	matrix ml_mst = (ml_mst \ tempmat)
	di "Iteration `i' complete."
}
mata : st_matrix("ml_mst_var", mm_colvar(st_matrix("ml_mst")))
matrix ml_mst_var = diag(ml_mst_var)
local colnames : colnames ml_ate_mst
mat colnames ml_mst_var=`colnames'
mat rownames ml_mst_var=`colnames'

bootstrapres ml_ate_mst ml_mst_var nobs
est store predictuse

coefplot (totaluse, mcolor(gs4) msymbol(Th) ciopts(lcolor(gs4))) (predictuse, mcolor(black) ciopts(lcolor(black))), ///
		vertical xtitle("Months From/Since Weatherization") ytitle("Energy Usage (MMBtu)") ///
		xline(12.5, lcolor(black)) graphregion(color(white)) bgcolor(white) ///
		xlabel(, labsize(small)) legend(order(2 "Realized" 4 "Predicted (ML)")) ///
		ylabel(, glcolor(gray%10))	
graph export "$maindir/Results/Figures/ML_trueVSpredict.pdf", replace


********** projected (engineering) energy savings
* Separating usage pre and post WAP, and creating artificial treatment variable
replace newsamp = 2 if newsamp==0
drop job_obs
sort Household newsamp
by Household : gen job_obs = _n
gen usage = ExistingConsumption/1000000 if job_obs==1 & newsamp==1
replace usage = ProjectedConsumption/1000000 if job_obs==2 & newsamp==1
replace usage = usage/12 /* transforming yearly to monthly */
winsor2 usage, cuts(0.5 99.5) trim replace
gen log_usage = log(usage)
drop ww_treat
gen ww_treat = 0 if job_obs==1 & newsamp==1
replace ww_treat = 1 if job_obs==2 & newsamp==1
label values ww_treat treat

* regression for projected energy savings - logs
reghdfe log_usage i.ww_treat , absorb(i.Household) vce(cluster Household)
nlcom exp(_b[1.ww_treat])-1, post
est sto projected_simple

* regression for projected energy savings - levels
reghdfe usage i.ww_treat , absorb(i.Household) vce(cluster Household)
est sto projected_simple_levels

	
******* realization rates
gen constant = 1
reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) residuals(res1)
est sto simple
reghdfe log_usage i.ww_treat , absorb(i.Household) residuals(res3)
est sto projected

reghdfe total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , absorb(i.Household#i.end_month i.month_year#i.CountyID) residuals(res1_levels)
est sto simple_levels
reg tau_ml constant if treated==1 & twoyears_prepost==1 & meterdays_monthly!=.  & newsamp==1 , nocons
est sto ml_levels
reghdfe usage i.ww_treat , absorb(i.Household) residuals(res3_levels)
est sto projected_levels

suest simple projected, vce(cluster Household)
nlcom (exp([simple]1.treated)-1)/(exp([projected]1.ww_treat)-1)
est sto att_realization
matrix rr1 = r(b)
scalar rr1 = rr1[1,1]
estadd scalar rr = rr1 : att_hhmonthcounty

est restore ml_simple
matrix beta = e(b)
sca mlbet = beta[1,1]
est restore projected_simple
matrix beta = e(b)
sca projbet = beta[1,1]
sca rr2 = mlbet/projbet
estadd scalar rr = rr2 : ml_simple

suest simple_levels projected_levels, vce(cluster Household)
nlcom [simple_levels]1.treated/[projected_levels]1.ww_treat
est sto att_realization
matrix rr1 = r(b)
scalar rr1 = rr1[1,1]
estadd scalar rr_lev = rr1 : att_hhmonthcounty_levels

suest ml_levels projected_levels, vce(cluster Household)
nlcom [ml_levels_mean]constant/[projected_levels]1.ww_treat, post
est sto ml_realization
matrix rr2 = r(b)
scalar rr2 = rr2[1,1]
estadd scalar rr_lev = rr2 : ml_simple_levels


esttab projected_simple ml_simple att_simple att_hhmonth att_hhmonthcounty ///
				using "$maindir/Results/Tables/ATT_fulltable.tex", ///
				rename(1.ww_treat a 1.treated a _nl_1 a constant a treated a) ///
				replace b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
				keep(a) coeflabels(a "IHWAP Treatment") nonotes ///				
				stats(rr N, fmt(%9.4g %9.0gc %9.4g) labels("Realization Rate" "Observations"))

esttab projected_simple_levels ml_simple_levels att_simple_levels att_hhmonth_levels att_hhmonthcounty_levels ///
				using "$maindir/Results/Tables/ATT_fulltable_levels.tex", ///
				rename(1.ww_treat a 1.treated a _nl_1 a constant a treated a) ///
				replace b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)  ///
				keep(a) coeflabels(a "IHWAP Treatment") nonotes ///
				stats(rr_lev N, fmt(%9.4g %9.0gc %9.4g) labels("Realization Rate" "Observations"))

